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ABSTRACT 

In a composite fluid system of two gravitationally coupled barotropic scale-free discs 
bearing a rotation curve v oc and a power-law surface mass density Eo oc r~ a 
with a = 1 + 2/3, we construct coplanar stationary aligned and spiral perturbation 
configurations in the two discs. Due to the mutual gravitational interaction, there are 
two independent classes of perturbation modes with surface mass density disturbances 
in the two coupled discs being either in-phase or out-of-phase. We derive analytical 
criteria for such perturbation modes to exist and show numerical examples. We com- 
pute the aligned and spiral perturbation modes systematically to explore the entire 
parameter regime. For the axisymmetric m — case with radial oscillations, there 
are two unstable regimes of ring-fragmentation and collapse corresponding to short 
and long radial wavelengths, respectively. Only within a certain range of the rota- 
tion parameter D 2 S (square of the effective Mach number for the stellar disc), can a 
composite disc system be stable against all axisymmetric perturbations. Compared 
with a single-disc system, the coupled two-disc system becomes less stable against 
such axisymmetric instabilities. Our investigation generalizes the previous work of 
Syer & Tremaine on the single-disc case and of Lou & Shen on two coupled singular 
isothermal discs (SIDs). Non-axisymmetric instabilities are briefly discussed. These 
stationary models for various large-scale patterns and morphologies may be useful in 
contexts of disc galaxies. 

Key words: waves — ISM: general — galaxies: kinematics and dynamics — galaxies: 
spiral — galaxies: structure — stars: formation. 



1 INTRODUCTION 

Rotating disc systems on various spatial scales are of 
broad astrophysical interest since most spiral galaxies, var- 
ious binary accretion systems, and proto-stellar and proto- 
planetary systems appear grossly in disc shape, which is 
believed to be a key intermediate stage that many astro- 
physical processes may attain (e.g., the phase of gas accre- 
tion onto a central black hole that drives an active galactic 
nucleus). It is therefore important to study the dynamical 
processes in disc systems for theoretical understanding and 
for astrophysical applications. Lin, Shu and co-workers pi- 
oneered the classic density wave theory in a differentially 
rotating disc system (Lin & Shu 1964, 1966; Lin 1987) and 
achieved a great deal of success in understanding the basic 
physics and dynamics of spiral galaxies (Binney & Tremaine 
1987; Bertin & Lin 1996) . The disc perturbation theory (lin- 
ear or even nonlinear) has proven to show broad potentials 
in dealing with problems of shapes and shaping, large-scale 



structures, instabilities in spiral galaxies and in other disc 
systems involving self-gravity, differential rotation and mag- 
netic fields (e.g., Binney & Tremaine 1987; Bertin & Lin 
1996; Balbus & Hawley 1998), as well as problems of angu- 
lar momentum transfer in accretion discs or planetary discs 
(e.g., Lynden-Bell & Kalnjas 1972; Goldreich & Tremaine 
1978). In many cases, perturbations developed in earlier 
stages prior to the moment when a system experiences more 
dramatic or violent processes (e.g., collapses), are crucial for 
dynamical evolution and are therefore worthwhile to pursue 
for their physical consequences. 

Among various problems in disc dynamics, a scale-free 
disc is often picked up by theorists for its relative simplicity 
and is explored as a powerful vehicle for a possible global 
analytical analysis. The term 'scale-free' here means that 
all physical quantities in the disc system vary as powers 
of cylindrical radius r (e.g., the linear velocity of disc ro- 
tation v oc r~P and the equilibrium surface mass density 
Eo oc r~ a ) with a and (3 being two related exponents. The 
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two examples in mind are the rigidly rotating Maclaurin 
discs and Kalnajs discs (Kalnajs 1972; Binney & Tremaine 
1987), where the angular rotation speed Q remains constant 
with v oc r. These discs are also known to have analytical 
normal mode spectrum, but are thought to rarely exist in 
nature. In contrast, thin discs with more or less flat rotation 
curves (i.e., v — constant) are common in most normal spi- 
ral galaxies as an important evidence for unseen masses of 
dark matter haloes associated with spiral galaxies. Besides 
these two limiting classes of discs with rigid and flat rota- 
tions, differentially rotating discs may have a rotation curve 
oc r _/3 with a rotation index (3 satisfying —1 < (3 < 1/2 and 
with (3 = 1/2 corresponding to a well-known Keplerian disc 
system.* 

Discs with complete flat rotation curves are usually re- 
ferred to as singular isothermal discs (SIDs), which form 
the simplest class in the family of scale-free discs. Since 
the introduction of SIDs by Mestel (1963), this idealized 
theoretical model has attracted a considerable attention in 
various astrophysical contexts of disc dynamics (e.g., Zang 
1976; Toomre 1977; Lemos, Kalnajs & Lynden-Bell 1991; 
Lynden-Bell & Lemos 1993; Syer & Tremaine 1996; Good- 
man & Evans 1999; Shu et al. 2000; Lou 2002; Lou & Fan 
2002; Lou & Shen 2003; Shen & Lou 2003; Lou & Zou 2004; 
Lou & Wu 2004). In the SID model, both the angular rota- 
tion speed fi and the equilibrium surface mass density Eo 
scale as r _1 , which is a scale-free condition. 

There has long been a paradox or controversy regarding 
stability analyses of scale-free discs because of the singular- 
ity as r — > 0. Starting from Zang (1976), who investigated 
a stellar SID numerically and argued that a scale-free disc 
can support no normal modes unless central cut-outs were 
introduced to remove the central singularity and to pre- 
scribe inner boundary conditions. Evans & Read (1998a, b) 
adopted Zang's approach to construct power-law discs with 
central cut-outs and examined numerically discrete grow- 
ing normal modes in an 'isothermal' stellar disc (i.e. with 
a constant velocity dispersion). In contrast, Lynden-Bell & 
Lemos (1993) claimed the presence of a continuum of un- 
stable normal modes for an unmodified SID. By specifying 
the phase of a postulated reflection of spiral waves from the 
origin r = 0, Goodman & Evans (1999) could define dis- 
crete normal modes for an unmodified gaseous SID. More 
recently, Shu et al. (2000) investigated spiral density wave 
transmission and reflection at the corotation circle and spec- 
ulated that the swing amplification process (Goldreich & 
Lynden-Bell 1965; Toomre 1981; Binney & Tremaine 1987; 
Fan & Lou 1997) across corotation allows a continuum of 
normal modes while proper 'boundary conditions' may se- 
lect from this continuum a discrete spectrum of unstable 
normal modes. 

Besides normal modes analyses, stationary perturbation 
configurations or zero-frequency neutral modes are empha- 



* Scale-free disc solutions do exist for /3 in the range of —1/4 < 
fi < 1/2 for warm discs according to our analysis. 



sized as marginal instability modes in scale-free discs (e.g., 
Lemos, Kalnajs & Lynden-Bell 1991; Syer & Tremaine 1996; 
Shu et al. 2000; Lou & Shen 2003; Shen & Lou 2003; Lou & 
Zou 2004). It is believed that axisymmetric instabilities set 
in through transitions of such neutral modes (Lynden-Bell 
& Ostriker 1967; Lemos, Kalnajs & Lynden-Bell 1991; Shu 
et al. 2000). By using properties of zero-frequency modes, 
Shu et al. (2000) further claimed that logarithmic spiral 
modes of stationary configurations also signal onsets of non- 
axisymmetric instabilities, a result compatible with the cri- 
terion of Goodman & Evans (1999) for instabilities in their 
normal mode treatment. Recently, Lou & Shen (2003) ex- 
tended results of Shu et al. (2000) in a gravitationally cou- 
pled composite disc system of one gaseous SID and one stel- 
lar SID in a two-fluid formalism. Stationary coplanar config- 
urations were readily constructed in such a composite SID 
system. 

The objective of this paper is to construct scale-free 
stationary configurations in a two-fluid stellar and gaseous 
disc system with a more general rotation curve v oc r^ 13 and 
equilibrium surface mass density profile Eo oc r~ a using a 
barotropic equation of state. The departure from the SID 
model will cause stationary configurations to vary signifi- 
cantly from those derived previously (Lou & Shen 2003) in 
some circumstances. 

This paper is organized in the following way. In Section 
2, we describe the theoretical formalism, obtain the back- 
ground rotational equilibrium and present linear coplanar 
perturbation equations. In Section 3, we discuss aligned and 
spiral disturbances in details and derive analytical criteria 
for stationary configurations. We summarize our results and 
give discussions in Section 4. Specific technical details are 
contained in the Appendices for the convenience of refer- 
ences. 



2 TWO-FLUID FORMALISM 

We adopt the two-fluid formalism sufficient for large-scale 
stationary aligned and unaligned coplanar disturbances 
(Kalnajs 1973) in a background rotational equilibrium with 
axisymmetry (Jog & Solomon 1984a, b; Elmegreen 1995; Jog 
1996; Lou & Shen 2003; Shen & Lou 2003). In this section, 
we provide the governing equations for a two-fluid compos- 
ite disc system, composed of a stellar disc and a gaseous 
disc presumed to be razor-thin. Given qualifications and as- 
sumptions, equilibrium properties of both barotropic discs 
characterized by rotation curves v oc r _f3 and power-law 
surface mass densities Eo oc r~ a with different proportional 
constants are described. We then derive linear coplanar per- 
turbation equations in such a composite disc system. 

It is important to note that the fluid formalism is well 
suited for large-scale dynamical behaviours in a gaseous disc 
but is an approximation to describe the large-scale dynam- 
ics in a stellar disc. The latter would be more appropriately 
modelled by the coupled collisionless Boltzmann equation 
(i.e., Vlasov equation) and Poisson equation (e.g. Julian & 
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Toomre 1966; Lin & Shu 1966; Zang 1976; Nicholson 1983; 
Binney & Tremaine 1987; Evans & Read 1998a, b), where 
an equilibrium distribution function (DF) is perturbed and 
the coupled Vlasov-Poisson equations are linearized to give 
the density wave dispersion relation. The introduction of 
an isotropic effective pressure term to mimic random stellar 
motions in a collisionless system is justified by the quali- 
tative agreement between a hydrodynamic formalism and a 
DF approach for treating collective particle dynamical be- 
haviours (e.g., Berman & Mark 1977). One illustrating anal- 
ogy between fluid and DF approaches is perhaps the den- 
sity wave dispersion relation in the WKBJ regime, namely, 
(mQ — uj) 2 = k 2 — 27rGSojfc| + k 2 v 2 in a gaseous disc and 
(mCl — uj) 2 — k 2 — 27rGS |fc|F in a stellar disc where F is 
the so-called reduction factor. In general, F is determined 
by the specific form of the DF, tends to reduce the gravity 
response and functions as an extra pressure term. The Q pa- 
rameter in the Toomre criterion for local axisymmetric sta- 
bility (Safronov 1960) also bears strikingly similar forms for 
both gaseous and stellar discs where for the latter the radial 
stellar velocity dispersion mimics the sound speed. This pro- 
vides an empirical rationale for treating the stellar disc by 
a simpler fluid approximation when dealing with the global 
axisymmetric stability for a composite disc system although 
the results may be quantitatively modified when we treat 
a stellar disc using the more exact (and more formidable) 
DF approach. The major deviation of a DF approach from 
the fluid formalism may occur in handling the corotation 
and Lindblad resonances. Therefore, for behaviours near the 
resonances, we need to rely on the DF approach for a stel- 
lar disc. In the present context of constructing large-scale 
stationary configurations, such resonances do not arise and 
the simpler two-fluid formalism for a composite disc system 
suffices. 



is due to the gravitational potential <j> through Poisson's 
integral 

-GT,(r',ip,t)r'dr' 



<))(r,6,t) = (tdip 



(4) 



r' 2 +r 2 - 2rr> cos(V> -6)] 1/2 ' 

where E = E s +E 9 is the total surface mass density. In a disc 
galaxy, the gravitational potential associated with a dark 
matter halo plays an important role. We shall take that into 
account later in a composite system of two coupled partial 
discs. * In equations (1) — (4), E 1 is the disc surface mass 
density, u 1 is the radial component of the fluid velocity, f is 
the z— component of the specific angular momentum about 
the z— axis, and II 1 is the two-dimensional effective pressure 
in the barotropic approximation, <j) is the total gravitational 
potential expressed in terms of Poisson's integral containing 
the total surface mass density E = E s + E 9 in a composite 
disc. Here, we assume that the stellar and gaseous disks 
interact primarily through the mutual gravitational coupling 
on large scales (Jog & Solomon 1984a,b; Bertin & Romeo 
1988; Romeo 1992; Elmegreen 1995; Jog 1996; Lou & Fan 
1998; Lou & Shen 2003; Shen & Lou 2003; Lou & Zou 2004). 

A barotropic equation of state assumes the relation 
between the pressure LI and the surface mass density E, 
namely, 

n = KTT , (5) 

where K > (i.e., warm discs) and n > are two constant 
coefficients and subscripts s (stellar disc) and g (gaseous 
disc) are implicit. This directly leads to the definition of 
sound speed a (in a stellar disc the velocity dispersion mim- 
ics the sound speed) by 
dU 



2 "J-J-0 r .„ n _l 

a = -— = nKE, 



(6) 



which gives a cx Eg™ with n — 1 for an isothermal sound 
speed. 



2.1 Basic Nonlinear Two-Fluid Equations 

We approximate both discs as razor-thin (i.e., infinitesimally 
thin) discs and use either superscript or subscript s and g 
to denote associations with the stellar and gaseous discs, 
respectively. The large-scale dynamic coupling between the 
two discs is due to the mutual gravitational interaction. For 
large-scale perturbations, we ignore diffusive effects such as 
viscosity, resistivity and thermal conduction, etc. Then the 
set of coplanar fluid equations for the stellar disc and the 
gaseous disc can be written out using the system of cylin- 
drical coordinates (r, 9, z) in the z = plane, such as 



r dr 



(rEV) + 



1 d 



du 1 i du 
+ u 



dt 



dr 



f du 1 



dt +u dr 
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fdf 
06 
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1 PIT _ 90 
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dr 



(1) 



(2) 



(3) 



where i — s, g denotes the stellar or gaseous disc here and 
throughout this paper. The coupling between the two discs 



2.2 Rotational Equilibrium of Axisymmetry 

It is straightforward to derive properties of the background 
rotational equilibrium of axisymmetry for the two gravita- 
tionally coupled discs, with physical variables denoted by a 
subscript 0. Let us first clarify several basic properties of a 
composite system of two scale-free discs. The background 
equilibrium surface mass densities of the two fluid discs Eg 
and Eq take the power-law form of tx r~ a with a common 
a exponent yet different proportional coefficients, while the 
disc rotation curves v s and v g take the power-law form of 
cx r"* 3 with a common /3 exponent yet different proportional 
coefficients. The special case of f3 = gives two flat rota- 
tion curves with v s / v g being allowed in general (Lou & 



t The construction of a composite system of two partial discs 
are described in Section 5. In our notation, the potential ratio 
T = 4>o/{4'0 + where ij>0 stands for the equilibrium background 
potential arising from the two discs and $ stands for that arising 
from an axisymmetric dark matter halo. Syer & Tremaine (1996) 
used dimcnsionlcss ratio / = 3>/</>o instead. 
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Shen 2003). For the background equilibrium, we also have 
uq = ug = 0, j'o = rv s and j'q = rv g . By imposing these con- 
ditions in equations (1) — (4) [particularly radial momentum 
eqn. (2)], we obtain 



v s + anK a (T,o) 



s\n-l 2 . r , msxn-l d(po 

v g + anK g (Eg) = r— . (7) 



To compute the gravitational potential <j)o arising from the 
equilibrium total surface mass density 

E = E^ + Eg = a s r~ a + (8) 

where oq an d °"o are two constant coefficients, we simply 
take equation (A5) in Appendix A of Syer & Tremaine 
(f996) and readily obtain 

4> = -27rGr(E^ + Eg)Po , (9) 
where we introduce an auxiliary parameter function 
r(-a/2 + l)r(a/2 - 1/2) 



(10) 



2r(-a/2 + 3/2)r(a/2) 
(Kalnajs 1971). 

The requirement of radial force balance (J7J for all radii 
(i.e., the scale- free condition) implies 

2/3 = a(n- 1) = a- 1 , (11) 



which gives the relationship among a, f3 and n, namely 

1+4/3 



a = 1 + 2/3 



and 



(12) 



1 + 2/3 

(Syer & Tremaine 1996). ft follows from n > in barotropic 
equation of state for warm discs that /3 > —1/4. For cold 
discs (i.e., K — > 0), this inequality is unnecessary. 

As discussed in Syer & Tremaine (1996), mass distri- 
butions with f3 > 1/2 (a > 2) would be unphysical be- 
cause they contain infinite point masses. Furthermore, Pois- 
son integral |g) for <^o converges for 1 < a < 2 and thus 
< (3 < 1/2 for a system of axisymmetry; the range of 
a (and thus of /3) is broader for nonaxisymmetric systems. 
However, the total force arising from axisymmetric equilib- 
rium surface mass densities remains finite in an extended 
range of /3 £ ( — 1/2,1/2) (0 < a < 2). In summary, we 
therefore have —1/4 < /3 < 1/2 [the left bound is implied 
by n > for warm discs and the right bound is required 
such that the central point mass will not diverge (Syer & 
Tremaine 1996)]. For cold discs (i.e., K = 0), the /3 range 
can be extended to -1/2 < /3 < 1/2. When (3 = for fiat 
rotation curves, we have surface mass densities proportional 
to r _1 corresponding to a composite system of two SIDs 
(Lou & Shen 2003; Shen & Lou 2003; Lou & Zou 2004; Lou 
& Wu 2004). 

According to equilibrium condition Q, we have 

V 2 S + A 2 = V 2 +A 2 g = 2nG(2(5Vo)r 1+2l3 (^ s + Eg) , (13) 

where v = Vr~ and a 2 = A 2 r' 20 /(1 + 2/3) with V and A 
being two constant coefficients. By introducing V = AD to 
define a dimensionless parameter D, we obtain 

A %Dl + l)r-< 1+2 « 



So 



ys 



2ttG(2/3Po)(1 + 5) ' 
A 2 g (D 2 g + l)Sr- ( -™ 
2ttG(2/3Po)(1 + 5) 



(14) 



where <5 = Eg/Eg is the ratio of the surface mass density of 
the gaseous disc to that of the stellar disc. We note that the 
value of 2(3Vo falls within (0, do) for /3 £ (-1/4, 1/2) and is 
equal to 1 when /3 = for the case of SIDs. 

An equivalent version of requirement (1131 is 



Al{D 2 s + 1) = A 2 g (D 2 g + 1) 



(15) 



where D s and D g are two dimensionless rotation parameters 
and rj = A 2 /A 2 = a 2 /a 2 is the square of the ratio of the ve- 
locity dispersion in the stellar disc to the sound speed in the 
gaseous disc. Note that A is actually related to the sound 
speed a tx r - * 3 [but scaled by a factor (l+2/3) 1//2 ] and the pa- 
rameter D is essentially the effective Mach number for disc 
rotation. We are going to express other equilibrium physical 
variables in terms of A and D. Besides, we have also intro- 
duced two dimensionless parameters to compare properties 
of the two discs. The first one is the surface mass density 
ratio 8 = Eq/Eq. The second parameter is the square of the 
ratio of the effective sound speeds in two discs r\ = A 2 /A g . 
For disc galaxies, ratio S can be either greater or less than 1 
depending on whether the system is stellar matter dominant 
or gas material dominant (in the early universe). Without 
loss of generality, we may take rj > 1 as the situation is sym- 
metric for rj < 1 and typically the stellar velocity dispersion 
(mimicked by a sound speed) in the stellar disc is greater 
than the sound speed in the gaseous disc. The special case 
of T) = 1 should give some familiar results of a single disc 
except for an additional mode due to gravitational coupling, 
as we have already learned from the simpler case of two cou- 
pled SIDs (Lou & Fan 1998; Lou & Shen 2003). 

The specific z— component angular momenta (j'q and 
j'q) and the sound speeds (a s and a g ) of the two discs in an 
equilibrium state simply read 

f = AMr 1 -? , (16) 



a 2 = nK^l)- 1 = A 2 /[(l + 2/3)r 2 "] . (17) 

Similarly, the disc angular rotation speed Q = jo/r 2 and the 
epicyclic frequency k = [(2fl/r)d(r 2 Q.) / 'dr] 1 ^ 2 are expressed 
in terms of two dimensionless parameters A and D as 

a = A l D l r- x - 13 , Ki = [2(1 - /3)] 1/2 a , (18) 

and therefore we have djo/dr = (1 — (3)v — rn 2 /(2Q,) that 
simplifies the linear perturbation equations displayed in the 
next subsection. 

For the convenience of comparison and cross referenc- 
ing, we note that our chosen notations for parameters have 
counterparts in those adopted by previous authors (Lemos, 
Kalnajs & Lynden-Bell 1991; Syer & Tremaine 1996). In 
Lemos et al. (1991), their notations a and v stand for 

A 2 2 A 2 D 2 



(1 + 4/3)r 2 < 3 



,.20 



(a 2 + v 2 y/ 2 [1 + (1 + 4/3)!) 2 ] 1 / 2 
here. While in Syer & Tremaine (1996), their notation w 
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stands for 



(l + 2(5)D 2 

for a full disc with their / = 0. These various adopted no- 
tations are relevant to the case of a single disc. All authors 
arrived at the same prescription for an axisymmetric back- 
ground in rotational equilibrium. 



2.3 Equations for Linear Coplanar Perturbations 

For small coplanar perturbations in a composite disc sys- 
tem denoted by subscript 1 associated with relevant physical 
variables, the perturbation equations can be readily derived 
by linearizing the basic equations (1) — (4), namely 



(19) 



djl r*4 
dt 2Vt 



■u\ + fli 



9 ( » B i • * 



for coplanar perturbations in the stellar disc and the gaseous 
disc, with the total gravitational potential perturbation 
given by 

-G(E| + T,{)r'dr' 



0i(r,0,t) = i>dij) 



(20) 

'o \r" + r* — 2rr' cos(i/j — 9)] LI A 

Assuming a Fourier component form of exp[i(cut — m0)] pe- 
riodic in time t and in azimuthal angle 9 for all perturbation 
variables with m > 0, we write for coplanar perturbations 
in the stellar and gaseous discs in the forms of 



Ej = fi 1 (r)exp[i(ait — mO)] , 
u\ = f/ I (r)exp[i(c 1 jt — mff)\ , 
j\ = J»exp[i(tjt - m0)] , 



(21) 



with the total gravitational potential perturbation in the 
form of 



(f>i — V(r)exp[i(cjt — mff)\ 



(22) 



within the disc plane at z = (we discriminate the imagi- 
nary unit i and sub- or superscript i for two discs). By sub- 
stituting expressions 1211 — 1221 into equations 1191 — 1201 . we 
readily derive for the stellar and gaseous discs 



1 d 



■ P 







\(lo — mf2j)/i* H — (rY? U l ) — imEp 

r dr i - 

i(w - mft)^ - 2SU— = ~ + V 

V (IV \ z_/ri 



i(w - mflijr + g-tr = + » 

and for the total gravitational potential perturbation 
-G{n" + fi 9 ) cos(mtp)r' dr' 



(23) 



V(r) = fdtp 



(24) 



'o (r' 2 + r 2 — 2rr' cos 1/)) 1 / 2 

We now use the last two equations in 123H to express U and 
J in terms of $ = a 2 /i/Eo + V for the stellar and gaseous 



discs 

U* = 



(w - rnfli) 2 - k? 

1 

r 



2i2i h (w — rail,) — 

r or 



(oi — raft) 2 



(oj — raft) - 



Kf d 
2ft cfr 



(25) 

Substitution of expressions j^H into the first equation of 
<2.'-!> leads to 

= (w — mft)// + 

(W — ?7lft) 2 — ft 2 

mEj 



r dr 

- 2ft h (w -raft) — * 

r dr 

. . 772 k~ d 



(26) 



r[(w - mft) 2 - k 2 ] 

for the stellar and gaseous discs. 

Based on equation 12611 . we construct stationary per- 
turbation solutions with u> = 0. With the axisymmetric 
background equilibrium conditions derived in Section 2.2, 
we rewrite equation l|2(il by setting uj = 0, namely 

8 1 

n + 

1 fm 2 +2/3 



D 2 (l + 2/3)(ra 2 -2 + 2/3) \ r 
(l + 2/3)(l + D 2 3 ) V 



dr 



dr 2 



(27) 



rfj, + 



2ttG{2[3Vo 
for the stellar disc and 



) 1 + 6 



= 



2/3)(m 2 -2 + 2/3) 
; , , (l + 2/3)(l + £% 



-2/3 



VS 



dr 



, 



dr 2 



(28) 



2ttG{2(3P ) 1 + S 
for the gaseous disc, respectively. The above two equations 
l|27l and 128(1 are to be solved together with Poisson's inte- 
gral 12411 . Note that equations 1271 and 128H are valid only 
for m / 0. In order to investigate the axisymmetric ra = 
case, we should take a different limiting procedure by first 
setting m — in equation 12611 before letting uo — > (Lou & 
Zou 2004). 



3 ALIGNED AND SPIRAL CASES 

We investigate in this section stationary density wave pat- 
terns in an inertial frame of reference using equations 1271 
and 128H coupled with Poisson integral 12411 . We distinguish 
two types of coplanar disturbances, that is, aligned and un- 
aligned (logarithmic spiral) coplanar perturbations. Aligned 
perturbation patterns correspond to distorted streamlines 
with the maximum and minimum radii at different radial 
locations lined up in the azimuth, while unaligned or spiral 
perturbations correspond to distorted streamlines with the 
maximum and minimum radii shifted systematically in az- 
imuth at different radial locations (Kalnajs 1973). Aligned 
perturbations relate to purely azimuthal propagations of 
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density waves (Lou 2002; Lou & Fan 2002) while the spi- 
ral perturbations relate to both azimuthal and radial prop- 
agations (Lou 2002; Lou & Fan 2002; Lou & Shen 2003; 
Lou & Zou 2004). Furthermore for aligned cases, we con- 
sider perturbations that carry the same density power-law 
dependence as that of the background equilibrium disc does. 
In contrast, we consider logarithmic spiral perturbations 
for nonaxisymmetric spiral cases (Kalnajs 1971; Syer & 
Tremaine 1996; Shu et al. 2000; Lou 2002; Lou & Fan 2002; 
Lou & Shen 2003; Lou & Zou 2004; Lou & Wu 2004). 



3.1 Aligned Perturbation Configurations 

The aligned m = case is somewhat trivial in the sense 
of a re-scaling of the axisymmetric background equilibrium 
state (Shu et al. 2000; Lou 2002; Lou & Shen 2003; Lou & 
Zou 2004). For m > 1, we consider aligned perturbations 
that carry the same radial power-law dependence as that 
of the axisymmetric background equilibrium disc does. The 
perturbed surface mass densities and the total gravitational 
potential read* 



V 



-(1+2/3) 



(j 9 r -(l+2/3) 



(29) 

-2KGr(n' + fi B )V m (P) , 
where a s , a 9 are small constant coefficients and the param- 
eter function V m ((3) is defined by 

r(m/2 - + l/2)r(m/2 + (3) 



V m ((3) 



(30) 



2r(m/2- ( S+l)r(m/2+/3 + l/2) ' 

with -m/2 < P < (m + l)/2 (Qian 1992; Syer & Tremaine 
1996). The prescribed ranges of /3 £ ( — 1/4,1/2) for warm 
discs and of (3 € ( — 1/2, 1/2) for cold discs happen to satisfy 
this requirement for m > 1. Note that for the isothermal 
case of ft — 0, we simply have V m = 1/m which § is just 
the case of Shu et al. (2000), Lou (2002), Lou & Fan (2002), 
Lou & Shen (2003), and Lou & Zou (2004). One can readily 
derive the recursion relation in m of V m ((3) for a fixed (3 
value, namely 

V m +i((3)V m ((3) = [(rn + 2(3)(m + l-2(3)]~ 1 . (31) 

In both ranges of (3 G (-1/4, 1/2) and (3 E (-1/2, 1/2), it is 
also useful to derive the asymptotic expression of V m ((3) 



V m ((3) « (m 2 + 2(3 - 4(3 2 ) 



-1/2 



(32) 



for m > 2 with an accuracy better than 2%. Larger values 
of m would lead to higher accuracies. 



By imposing condition 1291 with m > 1, equations 1271 
and can be cast into the following forms, namely 



M = 



m 2 + 2/3 d d 2 , .. 

2- r— - [Hirfi + Gir^ 

r ar dr z 



g . m +2(3 d d \(tt g , r> 
» = I " 2 d^~ r d^< iH2r ^ + G2r/i 



(33) 



where the four coefficients Hi, H2, Gi and G2 explicitly 
involved are defined by 



Hi 



H 2 



Gi 



G-2 



1 



D1{l + 2(3){m? -2 + 2(3) 

(1 + 2<3)(1 + D 2 ) 



2(3V 



(1 + 5) 



D g {\ + 2(3) (m 2 -2 + 2(3) 



(l + 2/3)(l + g 2 ) V m 8 



2(3V 



(l + <5) 



(34) 



1 + Di 



Vn 



Di(2/?P )(m 2 -2 + 2(3) {1 + 5) ' 



P m 8 



D|(2/3P )(m 2 -2 + 2(3) (1 + 5) ' 



where Vo appears due to background variables C1J. Note 
that Vo diverges at (3 = and (3 = 1/2, and approaches zero 
as (3 -> -1/2. Meanwhile, 2(3Vo = 1 is regular at (3 = 0. By 
expressions 129H . equations 1331 can be rearranged into 



[1 - H^m 2 + 4(3 - 4(3 2 )]^ s 
[1 - H 2 (m 2 + 4(3 - 4/3 2 )]/i 9 



Gi(m 2 + 4P-4(3 2 )ii 9 , 

9 1 (35) 

G 2 {m 2 +4(3-4(3 2 ) i i s . 



We introduce below several handy notations for parameters 
that depend only on m and (3 to greatly simplify analytical 
expressions, namely 

A m {(3) = m 2 + 4(3 - 40 2 , 
B m {(3) = (1 + 2(3) (m 2 -2 + 2(3) , 
C((3) = (1 + 2(3) /(2(3V ) , 
H m ((3) = C((3)V m ((3)A m ((3) + B m ((3) . 



(36) 



From these expressions, one immediately knows that pa- 
rameter C(0) decreases monotonically with increasing (3 and 
< C < 7r 2 /{4[r(3/4)] 4 } =* 1.094 for (3 G (-1/4,1/2) and 
C(0) = 1 by taking the limit of f3 -> 0. 

A straightforward combination of the two conditions 
in equation (I35H leads to the following stationary dispersion 
relation for nonaxisymmetric aligned coplanar perturbations 



t For aligned coplanar perturbations with a radial variation dif- 
ferent from that of the background equilibrium state, the per- 
turbation potential-density pair consistent with the Poisson in- 
tegral EHJ will be = crV-\ V = -2nGr(iM s + fl")Vm(X) 
where the superscript i = s,g for stellar and gaseous discs, re- 
spectively, numerical factor 7-' m (A) = F(m/2 — A/2 + l)T(m/2 + 
A/2 - l/2)/[2r(m/2 - A/2 + 3/2)r(m/2 + A/2)] and the A range 
of —rn + l<A<m + 2is required. Following the same procedure 
of analysis, we can construct a more broad class of stationary 
coplanar aligned perturbation solutions. 

§ We have assumed m > 0, otherwise \m\ should be used instead. 



(1 - HiA m )(l - H 2 A m ) = G1G2AL . (37) 

By substituting the expressions of Hi, H 2 , Gi and G 2 in 
equation 1341 into equation 1371 and using the background 
relation D 2 = i](D 2 + 1) — 1 where 77 = A 2 /A 2 , we obtain 
a quadratic equation in terms of y = D 2 for the stationary 
dispersion relation 137^ 

C 2 y 2 + Ciy + C Q = , (38) 
where coefficients C2, Ci and Co are functions of m, (3, 5 
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and 77, and are explicitly defined by 



C 2 ^BrnHmV 



Ci = 



Co 



(Am + Bm)(Um - Br, 



(1 + 5) 

(Am + B m )(Hm + B m S) 
(1+5) 
(Am + Bm)(H m - Bm) 



(39) 



+ (Am + Br, 



(1 



f S) 

+ Br, 



)(Hm+B m 8) 



(1+5) 

Given specified parameters of (3, 8 and rj, equation 13811 can 
be readily solved analytically for different values of m as 

-Ci ± (C 2 - 4C 2 C ) 1/2 
Vl ' 2 = 2C 2 ' 

where the determinant A = C\ — 4C%Co remains always 
positive for m > 1 as shown in Appendix A so that station- 
ary dispersion relation l|38ft has two distinct real solutions 
for D 2 . To illustrate the procedure, we choose several sets 
of parameters to numerically solve equation I38|l . 

The (3 = (or equivalently, n — 1) case has been thor- 
oughly studied as a composite system of two coupled SIDs 
(Lou & Shen 2003; Shen & Lou 2003; Lou & Zou 2004) where 
the surface mass density profiles scale as r~ 1 with flat ro- 
tation curves. In the present model, (3 is allowed to take 
on values in the range ( — 1/4, 1/2) for warm discs. For two 
representative examples, we choose (3 — —1/8 and (3 — 1/4 
corresponding to barotropic index n — 2/3 and n = 4/3, 
respectively. There are two more free parameters to be cho- 
sen: the first is the ratio of the surface mass density of the 
gaseous disc to that of the stellar disc 8. We choose three 
trial values of 8 = 1/4, 1 and 4 for three different cases. 
We simply assume 77 > 1 as the stellar disc has a relatively 
higher 'temperature' (i.e., higher velocity dispersion). We 
note when 77 — 1, the three coefficients C2, Ci and Co in 
equation 13811 turn out to be independent of 8. For 77 = 1 
and (3 = 0, the situation is reduced to two SIDs with the 
same sound speed (Lou & Shen 2003; Shen & Lou 2003). 

The representative solutions of equation l|38|l for differ- 
ent sets of parameters (3, 8, rj and m will be discussed more 
specifically later on. We here offer several remarks. Mathe- 
matically, two eigen-solutions for coplanar perturbations can 
be found due to the gravitational coupling between the two 
discs, but they may not always satisfy the physical require- 
ment D 2 > simultaneously. 

It is important to realize that the two branches of D 2 
solution of equation 138H . yi and 7/2, are both monotonic 
functions of rj (either monotonically increasing or monoton- 
ically decreasing) once other parameters are specified. This 
rule also applies to the spiral cases discussed later(see the 
proof in Appendix B). Therefore, once we know the solu- 
tions at rj = 1 and 77 — » 00 the entire solution ranges are 
qualitatively determined. More fortunately, the solutions at 
the two boundaries of rj have explicit analytical forms as 
shown below. 

For 77 = 1, we obtain two real solutions of quadratic 




Figure 1. Variations of coefficients Ai(j3), Bi(/3), CVi(/3), H\(0) 
as denned in equation 1361 and solutions of Y^\ m —\ and 
Yj 3 \m=l with 77 = 1 as obtained in equation I4UI in the (3 range 
of p e (-1/4, 1/2). 



equation 1381 in the forms of 
and 



\rA Am 

1 ~ B~ 



V B _ (1 - CVm)A„ 

*1 — 



Tin 



(40) 



This case is special since the second expression Yf is the 
result of one single disc and the first expression Yj 4 is extra 
due to the gravitational coupling. For the physical solution 
branch with D 2 > 0, Table 1 contains information of phase 
relationship between jiP and fj, s for the m = 1 case with 
different values of (3 £ (—1/4, 1/2). For Y± being physical, 
\x a and /x s are in-phase as a counterpart of a single disc, 
while for being physical, fj, 9 and /_t s are out-of-phase and 
this has no counterpart in the case of a single disc. 

Meanwhile in the limit of 77 — > 00, we solve quadratic 
equation 1381 to derive the following two solutions 

A _ Am[H m S + (1 - CVm)Bm] 



Yr 



B m Hm(l+S) 
(Am + B m )(H m 8 + Br, 
B,nHrn(l + 8) 



(41) 



and 
Y B = -1 

Apparently, the latter Y^ is unphysical for being negative. 

For the phase relationship between the two surface mass 
density perturbations fj, s and pP , we make use of equation 
I35H to derive 

JX» = 1 - H L A m = _ (D 2 Bm ~ Am)(l + 8) 

V s GlAm CVmAm(l + D3) 



: _ B m (l + 8) (Am + B m )(l + 8) 



(42) 



CT^mAm Cl^mAm(l ~\~ 

and then to calculate the phase relationship of surface mass 
density perturbations for stationary perturbation modes by 
inserting the value of D 2 solution obtained. 

3.1.1 The m = 1 Aligned Case 

The m = 1 case is somewhat special corresponding to 1 < 
CVi < 7r 4 /{12[r(3/4)] 8 } 1.596 for (3 g (-1/4, 0) and 
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(-l/4,-(2 1 /2 _i)/ 2 ) (-(2 1 / 2 - l)/2,0) (0,1/2) 



c 2 + 


+ 


1 -CPi 


+ 


Aitfi) 


+ + 






Wi(0) 


+ 


Yf\ m =l + 




Y B \ 1 

J X I'm, — 1 


+ + 




+ + 



Table 1. Signs for various parameters when /3 falls within the 
three contiguous subintervals of j3 £ (—1/4,1/2) for the m = 1 
case. The phase relationship between fj, 9 and fi s is only given 
for the physical solution branch with D 2 > and remains valid 
for all combinations of 5 and 77. Specifically, V 1 j4 | m=1 is positive 
for P e (-1/4, -(2 1 / 2 - l)/2), while Yf\ m=1 is positive for 
P e (-(2 1 / 2 - l)/2, 1/2). 



to < CVi < 1 for fi 6 (0, 1/2), respectively. We know 
that for a composite system of coupled two SIDs with fi = 
(i.e. flat rotation curves), the aligned m = 1 case imposes 
no restriction on the dimensionless rotation parameter D 2 
(Lou & Shen 2003). This situation changes qualitatively for 
/3 / 0. The additional freedom of fi parameter rules out 
that equation 13811 be automatically satisfied for arbitrary 
D 2 . We will see that to a certain extent, the aligned m — 1 
case is very similar to the spiral m = 1 case in form. The 
dependence of solutions y\ and 7/2 on the square ratio of 
sound speeds 77 and on the ratio of surface mass densities 8 
is distinctly different from those for the m > 2 cases. First 
we rewrite equations 1401 and I4H for m = 1 as 



Y A \ , - Al 

y a: 4/3(7*i$ + Bi 

r ao |m=l — — 1 + 



v b, _{l-CVi)Ai 



Hi 

Y B \ 1 - 

J 00 m — 1 - 



-1 , 



BiWi(l+<5) : 

where except for a removable singularity 1 ' at /3 = for 



=1, the coefficients _4i, Bi, CPi and Hi as well as 



the Dj solutions Yf 4 and Y\ as functions of /3 are displayed 
in Fig. 1 within the open interval fi G ( — 1/4,1/2). The 
sign variations of each parameters within the open inter- 
val fi € (—1/4, 1/2) are further summarized in Table 1 for 
reference. 

For the m = 1 case, the lower branch (either j/i or 2/2)" 
of the solutions to equation 1381 is always negative as can 
be seen later on and is therefore unphysical. It is possible 
for the upper branch of D 2 S solution to be positive for a 
specific range of 77, depending on the parameters /3 and 5 
(see the critical r] c discussed below). The variation within 
different parameter regimes is subtle for the special m = 1 
case as well as as for the phase relationship. And we only 

"I When /3 = 0, we haveWi = while Yf\ m—l given by equation 
141)1 remains finite because the numerator also vanishes. 
I When m = 1 , 3/1 is the upper branch and 3/2 is the lower branch 
for j3 < (C2 > 0), while the situation is reversed for fi > 
(C?2 < 0). However when m > 2, y\ always remains to be the 
upper branch and 1/2 always remains to be the lower branch as 
C 2 > within the range of all /3 G (-1/4, 1/2). 



need to consider the positive portion in the upper branch of 
D 2 S solution. 

As indicated in Fig. 1 and Table 1, we divide the open 
interval of— 1/4 < fi < 1/2 into three subintervals to analyze 
properties of aligned coplanar perturbations with m—l. 

★ Case I for -1/4 < fi < -(2 1/2 - l)/2 when t/i > is 
the upper branch 

As solutions yi and 7/2 of equation 138H are both mono- 
tonic functions of rj (see the proof in Appendix B), we can 
use the explicit j/i and 1/2 solutions at r\ — 1 [see solutions 
1H 1 and 77 — > 00 [see solutions 1411 1 to well bracket the 
ranges of D 2 value. For r\ = 1, we have 



2/i 



yi A | m= i 



Yi S U=i 



2/2 

while for r\ 



(\-CVx)Ax 



(43) 



<0 , 



2/i = Yi 



Hi 

00, we have 



= -1 



B1H1Q. + 6) 



> -1 



2/2 



V s I 

1 oc m—l 



(44) 



-1 



In this case, the lower 3/2 branch remains always negative 
and thus unphysical, while the upper y\ branch increases 
monotonically with increasing S and decreases monotoni- 
cally with increasing rj. For solutions I44II . there exists a 
critical r\ c beyond which y\ becomes negative for fixed val- 
ues of fi and 5. This r\ c can be explicitly determined by an 
analytical expression 

(l-CPiMi(l + <5) 

^ = ms + ii-cv^B, ■ (45) 

This expression 1451 remains valid only when rj c > 1 which 
further requires 



< S < 



(CTi - l)Bi 
Hi 



<CVi-l, 



(46) 



(47) 



that in turn defines a critical value 5 C = (CVi — l)Bi/7*i. 
For 5 > S c , the upper 3/1 branch remains always positive and 
physical. 

By equation 1421 , the phase relationship for surface mass 
density perturbations with m = 1 corresponding to the 
physical portion of the j/i branch is given by 

El = _i _ Bi(l + S) 4/3(1 + S) 
M s CV1A1 CViAi(l + D1) ' 

which decreases monotonically with increasing D 2 in the fi 
subinterval of case (I) and thus increases monotonically with 
increasing r\ for the y\ branch. We therefore have 

tl { 0- + S - CVi)/(CVi) < 0, ifO<<5<<5 c 
< 11 s < \ -CViAi5/(Hi5 + Bi) < 0, if 5 > S c 

(48) 

where the left-hand bound corresponds to 77 = 1 and the 
right-hand bound corresponds either to 77 = r\ c when con- 
dition 14fciH is met or to 77 — » 00 when 5 > 5 C . At any rate, 
the phase relationship of surface mass density perturbations 
in the two coupled discs for the upper 7/1 solution branch 
remains always out-of-phase. 
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(b) Case II: -(2 1 /2 _ l)/2 < /3 < 



k m=1 , p=-0.24, 8=4 


... .. . .... 


m=1, p=-0.24, 8=1 


' "\ 


m=1 , |3=-0.24, 8=1/4 

- 






2 3 4 5 


6 7 8 9 1 


(a) Case I: -1/4 < /3 < 


-(2 x /2 _ i)/ 2 . 




m=1, p=-1/8, 8=1/4 






m=1, p=-1/8, 8=1 




m=1, p=-1/8, 8=4 











*V~* m=1, (3=1/4, 8=1/4 




m=1, p=1/4, 8=1 




'~~"~----.__rn=1,P=1/4, 8=4 







(c) Case III: < P < 1/2 . 

Figure 2. Two solution branches of equation 1381 as functions 
of rj with combinations of m = 1, P = —0.24, —1/8, 1/4 and 
5 = 1/4, 1, 4. For the same set of parameters in each panel (a), 
(b), (c), each linetype corresponds to two solutions of D% for a 
range of 1 < rj < 10. 



* Case II for -(2 
upper branch 

For rj — 1, we have 



1/2 



l)/2 < P < when yi > is the 



2/i = ^i 



(l-CPi)^i 



5*1 [m=l 



J/2 

while for ?7 



Ai 



Hi 
< 



>0 , 



(49) 



Y A \ 

1 oo |m — 1 



Si 

oo, we have 

1 , 4/3(Hi<5 + Si) 



> -1 



SiHi(l + <5) • (50) 

2/2 = y£|m=i = -1 
In this case again, the lower 2/2 branch remains always nega- 
tive and thus unphysical. The upper y\ branch decreases 
monotonically with increasing either 5 or rj. The critical 
value of rj beyond which y\ becomes negative for fixed values 
of /3 and S is again determined by the same expression 14511 . 
This criterion remains valid only when rj c > 1 that further 
requires 

(CVi - l)Bi 



5 > 6 C = 



Hi 



>CVi-l 



(51) 



This then implies that for < S < 8 C , the upper 7/1 branch 
of D\ solution remains always physical for being positive. 

Similarly for the phase relationship between fi 9 and fj, s , 
we note that by expression 1471 . fj, 9 / 11 s increases monotoni- 
cally with increasing in the (3 subinterval of case (II) and 
thus decreases monotonically with increasing rj. We there- 
fore have 



5 > 



{l + 6-CVi)/(CVi) > , 
-CVx Ai5/(Hi5 + Bi) > 



if S > S c 
if < 5 < 6 C 

(52) 

where the left-hand bound corresponds to rj — 1 and the 
right-hand bound corresponds either to rj — rjc when condi- 
tion I51|l is satisfied or to rj — » 00 when < 5 < 8 C . At any 
rate, the phase relationship between ji 9 and ji s for the up- 
per 2/1 branch solution remains always in-phase with ji 9 / ' ji s 
smaller than 5. 

★ Case III for < f3 < 1/2 when 2/2 becomes the upper 
branch 

For rj = 1, we have 

{l-CVl)Ay 



2/1 



Y A \ , 

I\ \m — 1 



(53) 



2/2 = Y 1 \ m=1 = 
while for 77 



|m— 1 



Hi 

00, we have 

i/3(Hi8 



2/2 = Voo \m= 

or else, 

2/1 = Y^,\m= 
2/2 = K«,|m= 



if < 5 < -— - 



> 



■Bi) 



B1H1Q. + 6) 



< -1, 

B 
Hi 



, iN Si 
if > 



(54) 



4P(Hi6 + Bi) 
B1H1Q. + 6) 



Bi 
Hi 



In this case, the lower yi branch remains always negative 
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P = 


-0.24 


P = 


-1/8 


P = 


1/4 


m 


5 


5c 


Vc 




Vc 


Sc 


Vc 


1 


1/4 


0.401 


1.821 


0.522 


1 


0.809 


1 




1 


0.401 


/ 


0.522 


2.018 


0.809 


20.812 




4 


0.401 


/ 


0.522 


1.350 


0.809 


3.959 


2 


1/4 


/ 


2.514 


/ 


2.375 


/ 


2.033 




1 




1.813 




1.803 




1.809 




4 




1.556 




1.567 




1.664 


3 


1/4 




2.641 




2.278 




1.782 




1 




2.147 




1.947 




1.681 




1 




1.881 




1.752 




1.604 



while by solutions 1411 for rj — > oo, we have 



Table 2. The critical values of r\ c and 5 C for the three cases of 
j3 = -0.24, -1/8 and 1/4 with different values of m = 1, 2, 3 and 
5. A slash "/" means that the critical rj c or 8 C do not exist. Note 
that only for the m = 1 case does there exist such a critical <5 C 
that depends only on /3. 



and thus unphysical. The upper 2/2 branch decreases mono- 
tonically with increasing either 8 or rj. The critical value of 
rj beyond which y% becomes negative for fixed values of f3 
and 8 is also determined by expression 1451 . This criterion 
is valid only when r] c > 1 which further requires 
(CPi-l)Bx . 



8 > S c 



Hi 



(55) 



this inequality in turn implies that for < 8 < 8 C the upper 
U2 branch remains always physical for being positive. 

Similarly by equation 1471 for the phase relationship be- 
tween fi 9 and fi s , the ratio \l 9 I jf decreases monotonically 
with increasing D a in the f3 subinterval of case (III) and thus 
increases monotonically with increasing rj. We therefore have 



5< 



M» f (1 + S - CVi)/(CVi) > , if <5 > 5 C 
M a < I -(CViAi5)/{Hi6 + Si) > , if < S < 5 C 

(56) 

where the left-hand bound corresponds to n = 1 and the 
right-hand bound corresponds either to rj = r\ c when condi- 
tion 1551 is satisfied or to rj — > oo when < 8 < S c . At any 
rate, the phase relationship between jj? and /i s for the up- 
per i/2 branch solution remains always in- phase with ji 9 ' j ' ji s 
greater than 8. 

Relevant details of illustrating examples for three cases 
of = —0.24, f3 = -1/8 and /? = 1/4 respectively are 
shown in three panels (a), (b) and (c) of Fig. 2 and are also 
summarized in Table 2. 



2/1 



v a _ A m [n m 8+ (1 -CV m )B m ] 



B m H m (l + 8) ( 58 ) 
2/2 = = -1 . 

As the limiting situations for rj = 1 and r; — > oo well bracket 
possible ranges of j/i and j/2 branches (see Appendix B), it 
is obvious that the upper j/i branch remains always posi- 
tive and the value of y\ increases with increasing 8 and de- 
creases with increasing either m or r\. Meanwhile, the lower 
2/2 branch has a specific critical value r\ c of rj beyond which 
the ?/2 solution becomes unphysical for being negative. This 
critical value rj c is given by 

{1 - CP m )A m (l + 8) 



Vc = 1 + 



(59) 



H m S + (1 - CP m )B m 
where we have Jim always greater than m for all m > 2 
within f3 £ ( — 1/4, 1/2). Such a critical r\ c always exists for 
the j/2 branch for all values of 8. In other words, there does 
not exist a critical value S c for 8 as in the aligned m = 1 
case. 

This critical value j/ c of rj decreases with increasing 
8, much like the case of two coupled SIDs investigated re- 
cently (Lou & Shen 2003). For physical regimes of y\ and 2/2, 
they both decrease monotonically with increasing rj, while 
2/i branch increases monotonically and 2/2 branch decreases 
monotonically with increasing 8. 

For the phase relationship between ji 9 and ji s , it is 
straightforward to show that the ratio ji 9 / 'jj, s decreases 
monotonically with increasing D 2 a and thus increases mono- 
tonically with increasing rj for m > 2. For the upper 2/1 
branch, we obtain 

^ fJ> CT^rnAm8 

jj. s H m 8 + B 
where the left-hand bound corresponds to rj — 1 and the 
right-hand bound corresponds to rj — » 00. In the specified 
range of f3 6 ( — 1/4,1/2), the ratio fi 9 / ji s remains always 
negative for the upper 2/1 branch, indicating surface mass 
density perturbations in two discs are out-of-phase. 
For the lower 2/2 branch in parallel, we derive 



(60) 



5< ^ < cv„ 

where the left-hand bound corresponds to rj = 1 and the 
right-hand bound corresponds to the critical rj c which makes 
Ds =2/2 = 0. Apparently, the lower 2/2 branch (if physical) 
means surface mass density perturbations in the two coupled 
discs are in-phase. 



1 + 5 



(61) 



3.1.2 The m > 2 Aligned Cases 

When m > 2, we have C 2 > 0, < CVm < 1, Am > 0, 
B m > 0, Um > in the open (5 interval j3 G (-1/4, 1/2). 
Therefore in this case of m > 2, yi and 2/2 remain always 
upper and lower branches, respectively. By solutions 1401 for 
rj = 1, we have 

V A Am n 

yi = Y 1 = — > , 

t5m 



V2 = Y^ = ^^ >0, 

rim 



(57) 



3.2 Spiral Coplanar Perturbation Configurations 

Stationary surface mass density perturbations in both discs 
scale in the forms of oc fie~ lm0 in azimuthal angle 9. For 
aligned perturbations, we have further taken ji oc r~ e where 
e is a positive/negative constant exponent. For example in 
subsection 3.1, we have chosen e = a = 1 + 2/3 for copla- 
nar perturbations carrying the same radial power-law de- 
pendence of the background equilibrium disc system. On 
the other hand, for e being a complex constant exponent, 
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1 1 

(c) m = 3, /3 = -1/8, <5 = 4, 1, 1/4 (d) m = 3, /3 = 1/4, 5 = 4, 1, 1/4 

Figure 3. Aligned solution curves yi and 3/2 as functions of 77 for different azimuthal periodicities m = 2, 3 and surface mass density 
ratio 5 = 1/4, 1,4 with two fixed values of fi = —1/8 and /3 = 1/4. For the same set of parameters in each panel (a), (b), (c), (d), each 
linetype corresponds to two solutions of for a range of 1 < r\ < 10. 



perturbations would appear in spiral forms, namely, the so- 
called logarithmic spiral <x r~' x ^ exp[— iSJ(e) In r] where 
5i(e) and Ss(e) are the real and imaginary parts of e. To 
ensure the gravitational potential perturbation arising from 
this perturbed surface mass density as computed by Poisson 
integral (|1J being finite requires — m + 1 < K(e) < m + 2 
(Qian 1992). Without loss of generality, we assume a set 
of logarithmic spiral density perturbations and the result- 
ing gravitational potential perturbation in a mathematically 
consistent manner** (Kalnajs 1971; Syer & Tremaine; Shu 



et al. 2000; Lou 2002; Lou & Fan 2002; Lou & Shen 2003; 
Lou & Zou 2004; Lou & Wu 2004). Specifically, we write 

fi s = a 3 r~ 3/2 exp(i£ In r) , /1 9 = <r 9 r~ 3/2 exp(i£ In r) , 
\/ = -27rGr( M s + M 9 )A^(C) , 

(62) 

where a s and a 9 are small constant coefficients, 

K r (n = r(m/2 + g/2 + l/4)r(m/2 - ig/2 + 1/4) 
mlU 2r(m/2 + iC/2 + 3/4)r(m/2-iC/2 + 3/4) 1 ' 



** In parallel with the aligned case of coplanar perturbations, 
potential-density pair 1621 is not the only available potential- 
density pair that satisfies the Poisson integral 1201 - For coplanar 
logarithmic spiral perturbations, a more general class of allowed 
potential-density pairs that are consistent with the Poisson inte- 
gral OOJ is fjP = a3r- x cxp(i£ lnr), V = -2irGr(jJ, s +/t 9 )£m(f, A) 
where the superscript j = s and g, respectively, the numerical 



factor £ m (£, A) = r(m/2 - A/2 + if/2 + l)r(m/2 + A/2 - i£/2 - 
l/2)/[2r(m/2 - A/2 + if/2 + 3/2)r(m/2 + A/2 - if/2)] and the 
A range of — m + l<A<m-|-2is required. Following the same 
procedure of analysis, we can construct a more broad class of 
stationary coplanar perturbation solutions for logarithmic spiral 
configurations in a composite disc system. 
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is the Kalnajs function (Kalnajs 1971) and £ is a kind of 
radial 'wavenumber'. We refresh a few properties of Mm(£). 
First, A/"m(£) is an even function of £ so only £ > will be 
considered later. Secondly, Nm(£) decreases monotonically 
with increasing £ > 0. Thirdly, < Nm < 1 for m > 1 while 
No is positive and can be greater than 1 for a sufficiently 
small £. 

The choice of such form of perturbations is different 
from that of Syer & Tremaine (1996), whose spiral pertur- 
bations were taken to be oc r -1-2 ' 3 exp(im£ lnr) for m > 
(analysis before subsection 3.4 in their paper) in our nota- 
tions. For axisymmetric stability analysis in their subsec- 
tion 3.4, they adopted the same spiral perturbations in the 
form of 1621 (see also Lemos et al. 1991). We note that our 
background equilibria as well as the adopted form of loga- 
rithmic spiral perturbations are themselves scale-free, sepa- 
rately, whereas combinations of the background equilibrium 
and perturbations are not scale-free except for the special 
(5 = 1/4 case (see Lynden-Bell & Lemos 1993). 

Parallelling with Vm for the case of aligned perturba- 
tions, there are two useful formulae for N m ((.) for logarith- 
mic spiral perturbations. The first one is the recursion rela- 
tion 



2 i ,-2i-l 



N m +l(tWm(.0 = [{m + 1/2Y + e 



(64) 



and the second one is the asymptotic expression for N m (£) 

N m (0 « (m 2 +C 2 + l/4)" 1/2 (65) 

for m 2 + £ 2 3> 1. For rn > 2, this asymptotic expression l|65|l 
is accurate enough to compute values of N m (£)- 

Using potential-density set of 1621 for logarithmic spi- 
rals, we rearrange stationary coplanar perturbation equa- 
tions l|27|l and l|28[l into the following forms with m > 0. 



M = 



rn 2 + 2/3 



2-?- - r^r | (#177/ + Gir/i 3 ) , 



d_ 

dr 
d 



d*_ 

dr 2 
d 2 



,r = {I!L±2£_ 2 ^_ P ^_],v /j ,^_ f v J , / r 



(66) 



where the four relevant coefficients Hi, H2, Gi and G 2 are 
defined here explicitly by 



Hi = 



1 



Df(l + 2/3)(m 2 - 2 + 2/3) 

(l + 2/3)(l + £> 2 ) N m 



Hz 



1 - 
1 



2(3V 



(1 + 5) 



D 2 (l + 2/3)(m 2 - 2 + 2/3) 



(l + 2/3)(l + D 2 ) NmS 



Gi = - 



(i + D'i) 



2/3P (1 + 5) 

Nm 



(67) 



D 2 (2f3To)(m 2 -2 + 2/3) (1 + 5) ' 



G, 



(l + D 2 9 ) 



NmS 



Dj(2(3P )(m 2 -2 + 2/3) (1 + 5) 



Rewriting equations 1661 with expressions 16211 for fi s and 



2/3) M 9 



(68) 



H 9 , we immediately obtain 
[l-tfi(m 2 + C 2 + l/4 + 2/3)K 

= Gi(m 2 + C 2 + l/4- 
[1 - H 2 (m 2 + f + 1/4 + 2/3)]fi 3 

= G 2 (m 2 + £ 2 + l/4 + 2/3)// . 

As for the aligned case in subsection 3.1, we define some use- 
ful notations for parameter combinations that will simplify 
our following derivations, namely 

A4/3,C) = m 2 + e 2 + l/4 + 2/3 , 
Bm(P) = (1 + 2/3) (m 2 - 2 + 2/3) , 
C(/3) = (l + 2/3)/(2/3P ) , 

H m (A5) =CNmAm+Bm ■ 

With convenient notations (16911 . equation 1681 leads to the 
following stationary dispersion relation 



(69) 



(1 - HiA m )(l - H 2 A m ) = GiG 2 ^ 2 



(70) 



for coplanar logarithmic spiral perturbations in a composite 
disc system. 

Substituting expressions (16711 of Hi, Hi, Gi and G2 into 
stationary dispersion relation 1701 and using the background 
condition D 2 = rj(D 2 + 1) — 1, we obtain one quadratic 
equation in terms of y = D 2 , namely 

G22/ 2 + Ciy + C = 0, (71) 

where coefficients G2, Gi and Go are functions of parameters 
m, /3, 8, rj and £, and are defined by 

G 2 =B m 'Hm'n , 

[Am + B m )(Hm - B„ 



Gl 



Go = 



(1 + 5) 

_ (Am + Bm)(Hm + Bm.8) 

(1 + 5) 
(Am + Bm)(Hm - Bm) 



(72) 



+ (Am + B„ 



(1 + 5) 

(Am + Bm)(Hm + BmS) 



(1 + 5) 

Given specific values for m, f3, 8 and rj, we readily 
solve quadratic equation 1711 analytically for each 'radial 
wavenumber' £ for stationary logarithmic spirals. Again for 
a non-negative determinant as proven in Appendix A, there 
are two real D 2 solutions to equation 1711 . namely 

-Gi±(G 2 -4G 2 G ) 1/2 



2/1,2 



2G 2 



In addition to the well-studied isothermal /3 = case 
(Lou & Shen 2003; Lou & Zou 2004), we follow a similar 
procedure of analyzing the aligned case to construct loga- 
rithmic spiral configurations for the composite disc system 
in the range of /3 £ ( — 1/4, 1/2). For astrophysical relevance, 
we specifically choose /3 = —1/8 and /3 = 1/4 as illustrat- 
ing examples. Values of 8 are chosen as 1/4, 1 and 4 with 
77 > 1. For rj = 1, coefficients G 2 , Gi and Go as defined by 
expressions l|72|l are again independent of S. 

We now derive y — D 2 solutions below for the cases 
of rj — 1 and rj — > 00 separately for the same rationale as 
explained in the analysis of the aligned case. In terms of 
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coefficients 1691 , the two explicit D% solutions to stationary 
dispersion relation 1711 when r] = 1 are 



■ktA An 

1 ~ B~ 



and 



B _ (1 - CM m )Ar, 

i — 



H n 



(73) 



Since the ratio of sound speeds in the two discs are the same 
for rj = 1, the composite system of two coupled discs may 
be treated as one single disc to a certain extent. In fact, 
the second expression Y\ in equation 1731 is simply the 
result for the case of a single disc, while the first expression 
Y-^ in equation 1731 is additional due to the gravitational 
coupling between the two discs. Under some circumstances 
(e.g., m — and 1) when Y/ 1 remains always negative, we 
may practically regard the two-disc system as being identical 
with the case of a single disc for 77 = 1. 

The two explicit D 2 solutions to stationary dispersion 
relation l|710 in the limit of r\ — » 00 are 

A _ AmpimS + (1 - CMm)B m ] 



Y, 



and Y~ 



Meanwhile, the phase relationship for surface mass den- 
sity perturbations reads from l|68|l as 

(jj>_ = 1 - giAn = , _ (DsBm -A m )(l + S) 
M s GxAn Af m A m (l + D'i) ■ ( > 

As expected, all the expressions for the logarithmic spiral 
case can be obtained from those for the aligned case by sim- 
ply replacing V m with J\f m resulting from different potential- 
density pairs. This comes naturally from the perspective 
that both aligned and spiral configurations are coplanar den- 
sity waves propagating relative to the discs in either purely 
azimuthal directions or both radial and azimuthal directions 
(Lou 2002; Lou & Shen 2003). 

3.2.1 Marginal Stability of Axisymmetric Disturbances 

It is reminded that for the aligned case, axisymmetric m — 
perturbations merely represent a rescaling of the background 
equilibrium state of axisymmetry. For the m = case with 
radial oscillations, we should not start from equations 1271 
and 1281 . but instead^ should use equation 1261 by first 
setting m — with ui 7^ and then take the limit of uo — > 0. 
It is then straightforward to obtain 

[1 - H^ 2 + 1/4)K = G^ 2 + 1/4K , 
[l-i/2(C 2 + l/4)K=G 2 (C 2 + l/4)/i s , 

where coefficients Hi, H2, Gi and G2 are evaluated by sim- 
ply setting m = in definitions 1671 . It is clear that equation 
(O for the m = case with radial oscillations differs from 
equation 1681 for non-axisymmetric spiral cases. In defini- 
tions 1691 . we only need to replace Ao by A'o = £ 2 + 1/4 

tt It turns out that the outcome of Shu et al. (2000) with j3 = 
was all right in this regard because the two different limiting pro- 
cedures happen to give the same dispersion relation. Lou (2002) 
made a mistake in this regard because magnetic field breaks the 
degeneracy for the two different limiting procedures (see subsec- 
tion 3.2 and corrections in Appendix B of Lou & Zou 2004). 
M For the isothermal /3 = case, the two limiting procedures lead 
to the same result similar to the SID case of Shu et al. (2000). 



and redefine Ho = CNoA'o + Bo, and then use the explicit 
solutions derived for m > spiral cases in Section 3.2 by 
setting m = 0. The marginal stability curves thus obtained 
are displayed in Fig. 4—9. We show the complete y = D 2 
solution structure versus £ while only the portions above 
y = are physically sensible. 

To analyze the solution properties as parameters /3, 8, 
rj and 'radial wavenumber' £ vary for the axisymmetric m = 
case, we use asymptotic expression 1651 for Na(£,) and 
then use recursion relation Ki ll to derive an approximate 
analytical expression to compute the value of Ao(£), namely 



AAo(C) 



l F(l/4 + iC/2)r(l/4-ig/2) 
2 r(3/4 + i£/2)r(3/4 - i£/2) 

(9/4 + C 2 )(49/4 + C 2 ) 



(77) 



(1/4 + C 2 )(25/4 + £ 2 )(65/4 + e) 1/2 
with relative errors less than 0.5%. 

In the f3 interval of /3 6 ( — 1/4, 1/2), the quadratic co- 
efficient C2 of stationary dispersion relation 1711 vanishes 
when 



Ho = CAfoA'o + B = 



(78) 



which for a specified (5 determines the value of £ = £ c at 
which the value of D 2 will approach infinity. It is found 
that in the open interval of j3 G ( — 1/4,1/2), there exists 
only one such £ c where TLo vanishes as can be numerically 
computed^. Moreover for < £ < £ c , we have 7io < so 
that C2 > 0, while for £ > £ C) we have 7io > so that 
C2 < 0. This means that for wavenumber £ < f e , the upper 
branch is always j/i, while for wavenumber £ > £ c , the upper 
branch is always 3/2- 

First, we closely examine the /3 = —1/8 case. This 
critical £ c is numerically determined from equation l|781 as 
£ c = 1.217. We choose rj — 1 as a limiting case 1 ' 1 ' and the 
two explicit D 2 solutions from equation 1731 are 

A'o _ _16 f2 _ £ 

Bo ~~ 27 27 ' 2 (79) 



Y? 



(1 - CNo)Ao (1 - CAfoW + 1/4) 



Yl Ho CAA (C 2 + 1/4) -27/16 ' 

We stress here that Y-f and Y^ 3 do not correspond to y\ 
and 2/2 solutions, respectively, in a simple manner, because 
signs of coefficients vary with f. For instance in Fig. 4, we 
display y\ in solid line and j/2 in dashed line. Across £ c , solu- 
tion structures changes abruptly. For physically reasonable 

One can see this by directly comparing equations 1761 and 1681 . 
which are identical for m = if (3 = 0. 

S§ Ho increases monotonically with increasing £ and attains its 
minimum value at £ = with H.o_min 

< for all 6 (-1/4, 1/2), 
implying only one such £c where Ho = 0. See Appendix C for 
details. Only for one exceptional case of /3 ~ —0.130, this £ c 
happens not to be a divergent point for D 2 . See Appendix D for 
details. 

■ ■ For n = 1, we find the two-disc case is effectively identi- 
cal with the single-disc case as one solution branch (i.e., YA = 
A'o /Bo) remains always negative and is therefore unphysical. The 
other solution branch [i.e., K = (1 — CAfa)A 0/H0] is simply the 
same solution of the single-disc case (Syer & Tremaine 1996). 
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m=0, |3=-1/8, r,=1 

ring fragmentation 

£ =1.217 


collapse ^\ 







collapse 



m=0, p=1/4, r,=1 



\ ring fragmentation 



\ =3.159 



0.5 



2 3 4 5 



Figure 4. Two branches of solutions to equation 1711 for m = 0, 

fi = — 1/8 and r\ = 1. The divergent point is £ c = 1.217. In this 
case of r\ = 1, the value of 8 can be arbitrary. The y± branch is 
plotted in heavy solid line segments and the j/2 branch is plotted 
in heavy dashed line segments. 



Figure 6. Two branches of solutions to stationary dispersion 
relation 1711 with radial oscillations for m = 0, /3 = 1/4 and 
•q = 1. The divergent point is at £ c = 3.159. In this case of r) = 1, 
8 can be arbitrary. 



marginal stability curves, we only need the portions above 
y — 0. The two unstable regimes shown in Fig. 4 are the ring 
fragmentation regime where a composite disc system rotates 
too fast to be stable and the collapse regime where a compos- 
ite disc system rotates too slowly to be stable against large- 
scale Jeans collapse (Lemos et al. 1991; Syer & Tremaine 
1996; Shu et al. 2000; Lou 2002; Lou & Fan 2002; Lou & 
Shen 2003; Lou & Zou 2004). These marginal stability curves 
can also be derived from the time-dependent WKBJ analy- 
sis by imposing the scale-free disc conditions (Shen & Lou 
2003), with the more straightforward D s —criterion equiva- 
lent to the effective Q parameters presented by Elmegreen 
(1995) and Jog (1996). By varying the sound speed ratio 77 
and disc density ratio S, we obtain similar marginal stabil- 
ity curves in Fig. 5. The trends are qualitatively the same as 
the isothermal fi — case (Lou & Shen 2003; Shen & Lou 
2003). In other words, a composite disc system is less stable 
as compared with a single disc system for overall axisym- 
metric instabilities but becomes more difficult for large-scale 
collapses (Lou & Shen 2003; Shen & Lou 2003). 

Next, we consider the case of (3 — 1/4. The divergent 
point now becomes £ c = 3.159. For qualitative results, we 
again start from the special case of 77 = 1 with the two 
D 2 solutions of stationary dispersion relation l|71|l explicitly 
given by 

A'o 
Bo 

(l-CJVo).A'o 



= -4r/9 - 1/9 , 

(l-CA/o)(£ 2 + l/4) 



(80) 



Ho CA/"o(e + l/4)-9/4 • 

The corresponding marginal stability curves are displayed 
in Fig. 6. We further explored variations of the marginal 
stability curves for different sets of parameters in Fig. 7. 

There exists a critical /3 C above which the collapse 
regime disappears even for the 77 = 1 (single disc) case when 
the collapse region is largest. This critical f3 c = 0.436 is de- 
termined by the condition of zero collapsed regime for the 



m=0, P=0.45, r,=1 



* ring fragmentation 



% =11.346 




Figure 8. Two Dj solution branches to stationary dispersion 
relation 1711 with radial oscillations for m = 0, f3 = 0.45 and 
•q = 1. The divergent point is at £ c = 11.346. In this case of 
rj = 1, the value of 8 can be arbitrary. 



maximum of the lower-left branch, namely 
(1 - C/VoM'o 



Ho 



= 



(81) 



In order to see this clearly, we take (3 = 0.45 and obtain 
marginal stability curves for 77 = 1 as shown in Fig. 8 where 
no collapse regime appears. This result is consistent with 
that of Syer & Tremaine (1996) as can be seen from their 
fig. 2 for the marginal axisymmetric stable curve in terms 
of their w = 1/[(1 + 2(3)D'l] as noted earlier near the end of 
subsection 2.2 for notational correspondences. 

With the analysis technique developed by Shen & Lou 
(2003), we can perform time-dependent WKBJ analysis for 
a composite system of two coupled scale-free discs described 
in Section 2. For the above three cases with fi — —1/8, 1/4 
and 0.45, we display contours of frequency lu 2 in terms of 
the effective radial wavenumber f = K = \k\r and the ro- 
tation parameter D 2 in the same figure of exact stationary 
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(a) m = 0, /3 = -1/8, <5 = 1/4, r? = 5 



(b) m = 0, /8 = -1/8, 5 = 1, T] = 5 




(c) m = 0, /3 = -1/8, <5 = 4, »7 = 5 



(d) m = 0, /3 = -1/8, 5 = 4, r) = 10 



Figure 5. Two solution branches of stationary dispersion relation 1711 with radial oscillations for m = 0, /3 = —1/8, r\ = 5, 10 and 
5 = 1/4, 1, 4. In each panel (a), (b), (c), (d), the segments above y = constitute the marginal stability curves as indicated. 



perturbation configurations in Fig. 9 for the r) = 1 case cor- 
responding to the single disc case. The zero-frequency lines 
(i.e., marginal stability curves) for both precise and WKBJ 
approximation accord well with each other for large radial 
wavenumber when the WKBJ approximation is valid; for 
small radial wavenumber, the WKBJ approximation breaks 
down and the two regimes differ significantly as expected. 

In this context, we note the axisymmetric stability anal- 
ysis by Lemos et al. (1991) in a single-disc system. Lemos et 
al. (1991) derived the same axisymmetric background equi- 
librium state for a single disc. For perturbations, they im- 
posed adiabatic approximation with the adiabatic index 7 
greater than 1 and independent of the barotropic index n 
used for the equilibrium state, which is different from Syer 
& Tremaine (1996) and our present analysis. 

In the global analysis on axisymmetric m = instabili- 
ties with radial oscillations, stationary uj — wave patterns 
mark the onset of instabilities (Lemos et al. 1991; Syer & 



Tremaine 1996; Lou 2002; Lou & Fan 2002; Shu et al. 2000; 
Lou & Shen 2003; Shen & Lou 2003; Lou & Zou 2004) in a 
composite disc system. Apparently, there are two unstable 
regimes, namely, the long wavelength collapse regime and 
the short wavelength ring fragmentation regime. Therefore, 
the stability criterion for Dg falls in a range whose width 
increases with increasing (3. Both regimes of the collapse in- 
stability and the ring fragmentation instability are reduced 
for larger values of j3. As already noted, for f3 > 0.436, the 
collapse regime disappears completely. For sufficiently small 
values of /3 < —0.130, the stable range of Dg does not exist 
(see Appendix D for details). As remarked earlier, a com- 
posite system of two coupled discs is less stable than a single 
disc system. The introduction of an additional gaseous disc 
with larger 5 and r\ will reduce the overall stable range of Dj. , 
while tending to suppress the regime of collapse for large- 
scale instabilities (Lou & Shen 2003; Lou & Zou 2004). 
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(a) m = 0, P = 1/4, <5 = 1/4, r) = 5 
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Figure 7. Two D 2 solution branches of stationary dispersion relation 1711 with radial oscillations for m = 0, P = 1/4, r? = 5, 10 and 
5 = 1/4, 1, 4. In each panel (a), (b), (c), (d), the curve segments above y = constitute the marginal stability curves as indicated. 



3.2.2 The Logarithmic Spiral m = 1 Case 

The logarithmic spiral m — 1 case behaves qualitatively 
different from the aligned m — 1 case. The reason is simply 
that the additional radial wavenumber parameter £ will also 
alter values of coefficients in equation 11711 . 

In this case, we again use asymptotic expression 1651 to 
estimate M(£) and then use recursion relation l|64^ to derive 
an approximate analytical expression for M(£)> namely 



M 



1 T(3/4 + ig/2)r(3/4 - ig/2) 
2r(5/4 + i£/2)r(5/4-i£/2) 

(25/4 + £ 2 )(65/4 + e 2 ) 1/2 
(9/4 + £ 2 ) (49/4 + £ 2 ) 



(82) 



with relative error less than 0.5%. With Ai = ? 2 + 5/4 + 2/3 
and Bi = 4/3 2 — 1 according to definitions 1691 . we immedi- 
ately obtain 



Hi 



CM (£ 2 + 5/4 + 2/3) +4/? 2 - 1 , 



(83) 



which increases monotonically with increasing £ for fixed ft 
values and attains its minimum value"" at £ = as 

Wijnin = Wi|«=o >0 for p £ (-1/4, 1/2) . (84) 

This condition l)84ft means that no such £ exists to give Ti\ — 
0. Together with other relevant inequalities that Ai > 0, 
Si < 0, Ai +Bi > 0, < CM < 1, Hi > and C 2 < hold 
for all ^ > within the open interval of f3 E (-1/4, 1/2), 
we know for sure that and j/i are the upper and lower 
branch D 2 solutions, respectively. For t] = 1, we therefore 



III For P £ (—1/4, 1/2), the first order derivative of Hi with re- 
spect to £ is at £ = and positive for all £ > 0, while the 
second order derivative of H.i with respect to £ is always posi- 
tive at £ = 0. Hence £ = is a minimum for Hi- While we have 
obtained this result using the exact expression of Afi, the ap- 
proximate expression 1821 also works. See Appendix C for more 
details. 
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2/i = *i 



Ai 
Bi 



2/2 = fl | m= i = 



<0, 



(85) 



Hi 



> 



while in the limit of 77 



oo, we obtain 



(86) 



(.4 1 +g 1 )(H 1 <5 + /3 1 ) 
BiWi(H-J) 

2/2 = Voo|m=l = -1 

if inequality 5 > —Bi/Tii holds, or else 

2/1 = vj|m=i = -1 , 

„ _ (^l+fil)(ttl* + Bl) 

if inequalties < 8 < —Bi/Tii hold otherwise . 

It then follows that the lower 7/1 branch remains always neg- 
ative while the upper 7/2 branch becomes positive if 

(l-CMMi(l + <5) 



i< n < r, c _ 1 + + (1 _ cM)Bi] , 

which further requires the following inequality 
(CM - 1)23! 



8 > S c 



Hi 



(87) 
(88) 



(b) m = 0, = 1/4, r) = 1 



to be valid. 

This case is nearly identical with the case (III) in the 
aligned m = 1 case presented in subsection 3.1.1, except for 
the replacement of Vi by A/i and different definitions for 
Ai, Tii. The surface mass density perturbations in the two 
coupled discs are always in-phase for such configurations. As 
examples of illustration, we plot several cases with (3 — —1/8 
and 1/4, 8 = 1/4, 1 and 4 and 77 = 1 and 5 in Fig. 10 in terms 
of y = D 2 S versus £. Because the lower yi branch remains 
negative, we only show the upper y 2 branch in Fig. 10. 
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Figure 9. A comparison of stability boundary for exact global 
stationary perturbation configurations obtained here for m = 
with radial oscillations and of the zero-frequency boundary in the 
WKBJ analysis (Shcn & Lou 2003) for three cases: P = -1/8, 1/4 
and 0.45 with 77 = 1. 



3.2.3 Logarithmic Spiral Configurations with m > 2 

It turns out to be much simpler for m > 2 cases because 
when m > 2, we always have inequalities Am > 0, B m > 0, 
< CAf m < 1, Tim > and hence C2 > valid for all ranges 
of parameters under consideration. This greatly simplifies 
the analysis, as yi and 7/2 remain to be upper and lower 
branches, respectively. Meanwhile, we have Yi A > Yf > 
and YA > > Y* = -1. 



2/1 



For 77 = 1, we determine 
- Y a_A. 
1 ~~ ~B 

(1 - CMm)A 



>0 , 



(89) 



1 = 

while for 77 

Vi=Y oa = 



Tim 

— » 00, we have 

Am[HmS+ (1 



> 



- CM m )Br, 



> 



B m H m (l + S) ' - ' (90) 

2/2 = Y£ = -1 . 

Therefore, the upper yi branch remains always positive, 
while the lower 7/2 branch first remains positive for small 
77 and then becomes negative for 77 greater than a critical ?7 C 
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Figure 11. Two solution branches (given by the same linetype) of stationary dispersion relation 1711 for logarithmic spirals with 
m = 2, (3 = —1/8, 1/4, 8 = 1/4, 1, 4 and r) = 1, 5. For r\ = 1, these solutions are independent of 6. The lower branches are the 
counterparts of the single-disc case. Note the variation in ordering of two solution branches as r\ changes. 



given explicitly by 

, {1 - CAf m )A m (l + S) 



T), 



(91) 



[H m 5 + {1 - CAT m )B m ] ' 

This critical r\ c always exists for any given 8 because tj c 
remains always greater than 1 as dictated by equation 191H . 

Entirely similar to the aligned m > 2 cases, the phase 
relationship between surface mass densities fi 3 and fi s for 
the upper y\ branch of Dg solution is 

CAfmAmS 

'(H m 5 + B m ) ' 



-1< — < 



(92) 



where the left-hand side corresponds to 77 = 1 and the right- 
hand side corresponds to r\ — > 00. This branch always has 



out-of-phase relationship between surface mass densities 
and fx" in the range of j3 G (-1/4, 1/2). 

Meanwhile for the lower 1/2 branch, the phase relation- 
ship between surface mass densities fj, 9 and fj, s is determined 
by 



5< £- < 



(1 + 6) 



(93) 



where the left-hand side corresponds to r\ = 1 and the right- 
hand side corresponds to 77 = r\ c where D$ = y 2 = 0. 
This branch always has in-phase relationship between sur- 
face mass densities fi 9 and /-t s in the prescribed ft and i) 
ranges. 

For purpose of illustration, we present in Fig. 11 a few 
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solution examples in terms of y = D%) as functions of 'radial 
wavenumber' £ for specific parameters m = 2, (3 = — 1/8 and 
1/4, 5 = 1/4, 1 and 4 and r\ = 1 and 5. 



4 SUMMARY AND DISCUSSIONS 

The analysis presented in this paper is a generalization and 
extension of the previous work by Syer & Tremaine (1996), 
Shu et al. (2000), Lou & Shen (2003) and Shen & Lou (2003). 
We have constructed both aligned and logarithmic spiral, 
scale-free, coplanar stationary perturbation configurations 
in a composite system of two gravitationally coupled discs. 
While highly idealized, we have in mind, at least concep- 
tually, is a system of spiral galaxy consisting of one stel- 
lar disc and one gaseous disc, with a barotropic equation 
of state. This problem may then have relevance to distri- 
butions of stellar mass and gas materials in a disc galaxy. 
Qualitatively, the two branches of solutions derived in this 
paper suggest two possible coupled perturbation modes (not 
necessarily stationary; see Lou & Fan 1998b) where surface 
mass density perturbations in the stellar disc and in the 
gaseous disc exhibit either in-phase or out-of-phase corre- 
lations. These two distinctly different classes of perturba- 
tion modes are mathematically allowable, although there 
might be some kind of prevalence related to initial condi- 
tions or other uncertainties. For observational diagnostics of 
disc galaxies, one can obtain non-axisymmetric stellar struc- 
tures in the optical band (e.g. Rix & Zaritsky 1995) and de- 
rive HI maps for the gaseous disc component (e.g. Richter & 
Sancisi 1994). These two maps in different wave bands may 
be compared to see whether the stellar and HI gaseous arms 
are roughly coincident or apparently interlaced. The real sit- 
uation may be even more complicated than this. Active star 
formation processes in the optical arms will further consume 
HI gas and the places where HI gas clumps will trigger more 
star formation activities. Depending on the level of these 
interrelated processes, a phase shift between optical arms 
and HI arms may not be easily interpreted in terms of the 
out-of-phase perturbation modes. Perhaps, the most cogent 
evidence for out-of-phase density perturbations would be a 
lopsided disc galaxy where the gaseous and stellar disc com- 
ponents are comparable and the lopsidednesses for the two 
components are opposite. In a broader perspective, the ideal 
two-fluid approach adopted here may be applicable to other 
two-component disc system where the two components can 
be treated as ideal fluids with different temperatures, e.g., a 
composite disc system of stars and dusts or a composite disc 
system composed of young massive stars and relatively old 
stars, or even to composite disc systems with more compo- 
nents. These different 'hot' and 'cold' fluid disc components 
are coupled in the overall disc dynamics and contribute to 
various structures in multi-band observations. 

In order to satisfy the scale-free conditions in our 
model, both the stellar and gaseous discs in an axisym- 
metric equilibrium state have rotation curves v oc and 
surface mass densities £o oc r _1_2/3 with a barotropic in- 



dex n = (1 + 4/3) /(l + 2/3) for the parameter regime of 
/3 6 ( — 1/4,1/2). For both cases of aligned and logarith- 
mic spiral perturbations, we derive sensible values of D 2 to 
support such neutral or stationary density wave modes in 
an inertial frame of reference. There are two classes of sta- 
tionary density wave modes in a composite system of two 
coupled discs in general; this is in contrast to one class of 
stationary density wave modes in a single disc system. We 
now summarize our main results below. 

(i) Aligned Stationary Configurations 

For aligned configurations, we focus on the coplanar dis- 
turbances whose surface mass densities carry the same cylin- 
drical radial variations of the background equilibrium state. 

The aligned axisymmetric m = case represents merely a 
rescaling from one axisymmetric equilibrium to a neighbour- 
ing one (Shu et al. 2000; Lou 2002; Lou & Fan 2002; Lou & 
Shen 2003; Lou & Zou 2004), except that the rescaling here 
happens in both discs simultaneously. 

In contrast to the eccentric m = 1 case in a composite sys- 
tem of two gravitationally coupled SIDs with /3 = (Lou & 
Shen 2003), the aligned m = 1 configurations in two coupled 
scale-free discs with /3 are not trivial. Only one branch 
of D 2 solution is physically sensible. For j3 > — (2 1/2 — l)/2, 
this branch of D 2 solution stands as the counterpart of the 
single-disc case and the surface mass density perturbations 
in the two discs are in-phase. For f3 < — (2 1 / 2 — l)/2, this 
branch of D 2 solution has no counterpart of the single-disc 
case and surface mass density perturbations in the two discs 
are out-of-phase. There may or may not exist a critical n c , 
determined by expression l|45|l . beyond which the D 2 so- 
lution becomes negative and thus unphysical, depending on 
the value of S. Specific classifications and analyses have been 
presented in subsection 3.1.1. 

For m > 2 cases, we have derived two branches of so- 
lution for possible values of the rotation speed parameter 
D 2 such that aligned stationary perturbation configurations 
are sustained in a composite disc system. Of the two D 2 
branches, one is always the upper branch and thus physi- 
cal for being positive and the coplanar surface mass density 
perturbations in two discs are out-of-phase (see Lou & Fan 
1998). This branch of D 2 solution has no counterpart in the 
case of a single disc. Meanwhile, the other branch of D 2 so- 
lution stands as the counterpart of the single-disc case with 
coplanar surface mass density perturbations in the two discs 
being in-phase. Furthermore, this second branch of D 2 solu- 
tion decreases with increasing r\ and may become negative 
and thus unphysical when r\ exceeds a critical value n c that 
varies with f3, m and S as seen from definition 159H . In con- 
trast to the aligned m — 1 case, this critical n c always exists 
for any values of 8. Details and specific examples can be 
found in subsection 3.1.2. 

(ii) Stationary Configurations of Logarithmic Spirals 
For unaligned or spiral coplanar perturbations, we con- 
sider global logarithmic spiral configurations (Kalnajs 1971; 
Syer & Tremaine 1996; Shu et al. 2000; Lou 2002; Lou & 
Fan 2002; Lou & Shen 2003; Lou & Zou 2004). 
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For the axisymmetric m = case with radial oscillations, 
we have determined the marginal stability curves of D 2 ver- 
sus the 'radial wavenumber' £ for various values of (5. The 
limiting case of n = 1 reduces to the case of single disc case 
as if the secondary mode due to the gravitationally coupling 
between the two discs were absent. The axisymmetric sta- 
bility criterion is expressed in terms of the stellar rotation 
parameter D 2 . Those systems that rotate too slowly will 
succumb to large-scale instabilities in the collapse regime, 
while those systems that rotate too fast will fall into the 
ring-fragmentation regime for short-wavelength instabilities 
(Safronov 1960; Toomre 1964; Shu et al. 2000; Lou 2002; 
Lou & Fan 1998, 2002; Lou & Shen 2003; Shen & Lou 2003; 
Lou & Zou 2004). The stable range of D 2 against overall 
axisymmetric instabilities is expanded for larger /3 values. 
When (3 > f3 c ~ 0.436, the large-scale collapse regime will 
disappear completely. On the other hand when f3 < —0.130, 
the system cannot be stable at all. The D 2 criterion for 
axisymmetric instabilities with radial oscillations presented 
here is entirely equivalent to the w parameter used by Syer 
& Tremaine (1996) for the the case of a full single disc. The 
composite disc system becomes less stable than the single- 
disc system. The overall stable D 2 range will diminish for 
larger 8 or n, while the large-scale collapse instability by 
itself tends to be suppressed. Specific examples and some 
analysis techniques are presented in subsection 3.2.1 and 
Appendixes C and D. 

For the logarithmic spiral case of m = 1 , the lower branch 
of D 2 solution is always unphysical for being negative and 
the limiting case of rj = 1 corresponds to just the single- 
disc case. For a given 'radial wavenumber' £, there may or 
may not exist a critical n c beyond which the D 2 solution 
becomes negative, depending on the value of 8. This case 
is almost identical with the aligned case (III) described in 
subsection 3.1.1, except for a substitution of V\ with M\ and 
redefinitions for A\ and Tii. Details and specific examples 
can be found in subsection 3.2.2. 

For logarithmic spiral cases with m > 2 , we have ob- 
tained analytical results almost in the same forms of the 
aligned cases. There are two possible values of rotation pa- 
rameter D 2 that can sustain stationary logarithmic spiral 
configurations in a composite system. Of these two pertur- 
bation modes, one has no counterpart in the single-disc case 
(Lou & Fan 1998) with the surface mass density perturba- 
tions in the two discs being out-of-phase, while the other 
is the counterpart of the single-disc case and has in-phase 
surface mass density perturbations in the two discs. The 
out-of-phase mode always exists while the in-phase mode 
disappears when r\ exceeds a certain critical value rjc that 
depends on (5, m, 8 and £ by equation (I9H . Specific exam- 
ples can be found in subsection 3.2.3. 

(iii) N on- Axisymmetric Instabilities 

For axisymmetric instabilities it is well established that 
neutral modes mark the onset of instabilities. We have ana- 
lytically constructed non-axisymmetric neutral (i.e. station- 
ary) modes for both aligned and logarithmic spiral cases. 
Do these neutral modes or stationary configurations mark 



the non-axisymmetric spiral instabilities? Based on trans- 
missions and over-reflections of leading and/or trailing spiral 
density waves across corotation in a time-dependent anal- 
ysis, Shu et al. (2000) speculated that the condition for 
stationary logarithmic spiral configurations in a SID would 
determine whether spiral density waves could be swing- 
amplified (Goldreich & Lynden-Bell 1965; Fan & Lou 1997). 
The criterion of Shu et al. (2000) for swing amplification is 
consistent with that of Goodman & Evans (1999) in their 
normal-mode analysis. It is appears worthwhile to pursue 
the criterion for the onset of non-axisymmetric instabilities 
in terms of a normal-mode analysis for a composite system 
that will be performed in a separate paper. 

(iv) Partial Discs 

All computations and discussions above deal with full 
discs. One can also impose an axisymmetric gravitational 
potential associated with a background dark matter halo 
of axisymmetry and ignore disturbances in the dark matter 
halo caused by coplanar surface mass density perturbations 
in the composite disc system,* * * that is, the composite sys- 
tem is composed of two partial discs (Syer & Tremaine 1996; 
Shu et al. 2000; Lou 2002; Lou & Fan 2002; Lou & Shen 2003; 
Shen & Lou 2003; Lou & Zou 2004). By defining a dimen- 
sionless factor < J- < 1 for the ratio of the gravitational 
potential arising from the two discs together to that of the 
whole system including the dark matter halo, one may fol- 
low the same procedure for full discs to analyze the problem 
of a composite system of two coupled partial discs. Practi- 
cally, what one needs to do is to replace all V m and N m with 
J'Vm and J-Nm, respectively. The dynamical effect of this 
background dark matter halo tends to suppress axisymmet- 
ric instabilities (Syer & Tremaine 1996; Shu et al. 2000; Lou 
2002; Lou & Fan 2002; Lou & Shen 2003; Shen & Lou 2003; 
Lou & Zou 2004). 

(v) Magnetized Discs 

By synchrotron radio observations, one can infer spiral 
magnetic field structures in nearby spiral galaxies (e.g., Lou 
& Fan 2003 and references therein) . It is believed that this is 
generically true for distant spiral galaxies as well. The pres- 
ence of magnetic fields in spiral galaxies should affect global 
star formation rates and thus influence the evolution of disc 
galaxies. Technically, the interesting problem of including 
magnetic fields can become quite involved. More specifically, 
a stellar disc is gravitationally coupled with a magnetized 
gaseous disc. There are two relatively simple geometries to 
model configurations of magnetic fields. The first one is the 
so-called isopedically magnetized configurations (Shu & Li 
1997; Shu et al. 2000; Lou & Wu 2004). The second one 
is to consider coplanar azimuthal magnetic fields with their 
strengths scaled as powers of cylindrical radius r (Lou 2002; 
Lou & Fan 2002; Lou & Zou 2004). Now with a more general 
scale-free disc system investigated in this paper, it would be 

* * * This simplifying approximation is crudely justifiable for 
high velocity dispersions in a dark matter halo. By numerical sim- 
ulations, velocity dispersions of dark matter halo are presumably 
of the order of a few hundred kilometers per second. 
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very interesting to model magnetic fields for both isopedic 
and azimuthal configurations in a more general composite 
disc system, a problem also to be presented in a separate 
paper. 
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APPENDIX A: TWO REAL D% SOLUTIONS 

We here show that equation 13811 always has two different 
real solutions. The determinant of equation 1381 reads 



A = Gl - iC 2 G = ( Am + B -° 
\ l + o 

where 



P = C2V 



ClT] + Co = c 2 7? 



Ci \ 2 _ cj - 4c c 2 
2c 2 J 4c 2 



(Al) 



(A2) 



and 

C2 = (B m + H m 3) 2 , 

Cl = 25{H m - B m f ~ 2B m H m {l + S) 2 , (A3) 
CO = (H m + B m 8f . 
It follows further that 

Ai = c? - 4c c2 = -im m H m (H m - B m ) 2 5{l + S) 2 . (A4) 

Now in the open interval of /3 € ( — 1/4, 1/2) and for m > 2, 
we have C2 > 0, Ai < and thus p > for all 77 > 1. We 
have thus proven that A > 0. While the proof procedure is 
slightly different for the m = 1 case, we also can show that 
A > 0. It follows that equation 1381 always has two different 
real solutions for D 2 corresponding to the upper and lower 
branches, respectively. 

The same proof procedure can be repeated for the log- 
arithmic spiral cases corresponding to stationary dispersion 
relation 171H . which always has A > for m > 1. Moreover 
in the same spirit, we can show A > for the axisymmetric 
m — case. The equal sign corresponds to (A'o + Bo) — 0. 



APPENDIX B: MONOTONIC FUNCTIONS 

We here provide proofs that the two D 2 solutions of dis- 
persion relation Ij38^ for the aligned case and of dispersion 
relation 1710 for logarithmic spiral cases are monotonic func- 
tions of ?7. Therefore, one can make use of the explicit solu- 
tions at 77 = 1 and in the limit of rj — > 00 to well bracket the 
D 2 solution ranges. 
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First, we rewrite the coefficients C 2 , Ci and Co as de- 
fined by either expressions 1391 for aligned perturbations or 
expressions (17211 for logarithmic spiral perturbations, namely 

C 2 = a 2 r) , 

C 1 = am + h, (Bf) 
Co = a r/ + b , 

where coefficients a 2 , 01, 61, ao and &o are determined by 
directly comparing expressions 1)390 or l)72|l of the actual 
coefficients C 2 , Ci and Co that appear in the main text. 

With the two D 2 solutions explicitly given by 7/1,2 = 
(-Ci ± A 1/2 )/2C 2 , we have 

d yi , 2 A 1/2 (CiC 2 - C 2 C{) + (C 2 A' - AC 2 ) 



dr) 



4C 2 A 1 /2 



(B2) 



where / denotes the derivative with respect to r\ and A = 
G\ — 4C 2 Co is the determinant that has been proven non- 
negative in Appendix A. 

We next consider to solve equation y[ 2 = dy\. 2 jdr\ = 
that is equivalent to the following condition 

477 2 a 2 Ooa 2 - aiMo + a b'{) = . (B3) 

By straightforward substitutions of the expressions of a 2 , ai, 
61, ao, bo, this condition turns out to be 

4 V 2 B m H m {Am + B,l ;f (n Z ~ Bm)H = , (B4) 

(1 + oy 

which gives only one solution 77 = for any given sets of 
{m, (3, 8} for aligned perturbations or {m, (3, 5, £} for log- 
arithmic spiral perturbations. Since dyi :2 /drj axe continuous 
functions of 77, it then follows that for 77 > 1, dy 1>2 /dr/ remain 
either always positive or always negative for any specified 
sets of {m, (3, 8} for aligned perturbations or {m, ft, 8, £} 
for logarithmic spiral perturbations. In summary, the two 
D 2 solutions t/i j2 must be monotonic functions of 77 once 
{m, (3, 8} for aligned perturbations or {m, (3, 8, £} for 
logarithmic spiral perturbations are specified. 



APPENDIX C: PROPERTIES OF Hm 

We here study the variation of 7i m with respect to £ in the 
logarithmic spiral case. For m > 0, we consider 

An = m 2 + £ 2 + 1/4 + 2(3 , 

B m = (1 + 2(3) (m 2 -2 + 2/3) , 

C = (1 + 2(3)/{2pV Q ) , 

7~tm — CA/m-4-m "f" ■ 

It is then straightforward to show 
dTCm 

~w 

where 



= CN m {y m A m + 2£) 



(CI) 



(C2) 



*m(6 = -[^(m/2 + ^/2 + l/4)+V(m/2-^/2 + 3/4) 



■ 7>(m/2 - i£/2 + 1/4) - V(m/2 + i£/2 + 3/4)] 



where 7/) is the digamma function (or 7/)— function) defined 
by tp(x) = d[lnF(x)]/dx with T(x) being the T— function. In 
the above derivation, we have used the following relation 

dAf m 



d* 



(C4) 



It is also straightforward to show that \I/m(£ > 0) is 
negative, ^I' m (0) = and therefore 

dUr, 



= 



5=0 



Furthermore, we have 
d 2 H 



d£ 2 
dTLm 



> for (3 G (-1/4, 1/2) 



5=0 



di 



(C5) 



(C6) 



(C8) 



> for (3 6 (-1/4,1/2) . 

Therefore, H m increases monotonically with increasing £ 
and attains its minimum value at £ = 0, that is, 

Hm.min = H m \i=0 = CJ\f m (0)Am(0) + B m . (C7) 

For the m = case with radial oscillations, we should 
make the following two replacements 

Ao -> A' = e 2 + 1/4 , 

Ho CNoA'o + Bo , 

and repeat the same procedure for m > cases. The pre- 
ceding results remain valid. In summary, for m > 0, TL m 
increases monotonically with increasing £ and attains its 
minimum value at £ = 0. More specifically, we have 

Ho.rn.in = Wo|e=o = GA/b(0)(l/4) + (1 + 2/3)(-2 + 2(3) < , 
Hi_ min = Wi| e =o = CM (0) (5/4 + 2/3) + 4/3 2 - 1 > , 
H = Hm|j=o = CA/m(0)^4 m (0) + S m > for m > 2 . 

(C9) 



APPENDIX D: 

We here discuss more specifically the m = case with radial 
oscillations, for which we have 

Ao = £ 2 + 1/4 > , 

So = (l + 2/3)(-2 + 2/3) <0 , 

C = (1 + 2/3)/(2/3P ) , 

Ho = CAM'o + Bo . 



(Dl) 



For r\ = 1 (equivalent to the case of a single disc) , the phys- 
ical D 2 solution is 



(D2) 



(C3) 



B _ (l-CMo)Ao . 
Yl Ho >0 ' 

the other D 2 solution is negative and thus unphysical. It is 
noted that for m > 1, we have < CN m < 1 for all £ > 
and (3 £ ( — 1/4, 1/2). In contrast, the situation of m — is 
different, that is, CNo can be either greater or smaller than 
1 and this complicates the analysis. 

According to Appendix C, we infer 

= CA/o^o < for £>0, (D3) 

at, 

and is equal to zero at £ = 0. This implies that CNo decreases 
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Figure Dl. Two branches of D 2 solutions to stationary disper- 
sion relation 1711 for m = 0, = —0.2 and r} = 1. The divergent 
point is £ c = 0.807. In this case of T] = 1, the value of <5 can 
be arbitrary. The system is inevitable to bear ring-fragmentation 
instabilities and therefore there is no stable range of D 2 . We also 
show the deviation from contours of local WKBJ analysis, which 
are consistent with the global analysis at large 'radial wavenum- 
ber' §. 



monotonically with increasing £ and reaches the maximum 
value at £ = 0, namely 

(CA/i)«~ = CA/b(0)| <x . f/3> ^ 0436 ; (D4) 

Therefore for /3 < /3 C ~ 0.436, there is a critical £ c / at which 
1 — CTVo vanishes. If this £ c / coincides with £ c at which 7io 
vanishes, then there is no divergent point for Yf 8 . For this 
reason, we now check the possibility of f c / = f e . 

If £ c / = £ c happens for a specific /3 co , the following two 
equations must be satisfied simultaneously, 

(CM>.A'o + Bb)|e e ,=« e = 0, 

which gives the relation of £ c in terms of /3 co , 

.A'o + Bo = d + 1/4 + (1 + 2/3) (-2 + 2/3) = (D6) 

and thus 

C c = (7/4 + 2/3-4/3 2 ) 1/2 . (D7) 

Inserting solution (113711 into the first equation of (D5), we 
obtain numerically 

f3 co ~ -0.130, Cc' = £c ~ 1.193 . (D8) 

From the above analysis, we realize that for /3 co < P < 1/2, 
the axisymmetric marginal stability curves are typical as dis- 
cussed in subsection 3.2.1. In contrast, for —1/4 < (5 < (3 co , 
the axisymmetric marginal stability curves are different and 
there is no stable range of D 2 against overall axisymmet- 
ric instabilities, a result also implied in Syer & Tremaine 
(1996) (see their fig. 2). To see this more clearly, we con- 
sider a specific case of j3 — —0.2. The divergent point of Yf 8 
is located at £ c = 0.807. The stationary configuration for 
r\ = 1 (equivalent to the case of a single disc) is shown in 
Fig. Dl. 
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